Measurement exploratory data analysis

DEVICES = [
    'beaglebone-fan',
    'beaglebone-compressor',
    'beaglebone-pump',
    'beaglebone-refrigerator',
    'mafaulda-a',
    'mafaulda-b'
]

DEVICE = DEVICES[1]
T_WAVEFORM = 5  # (1 = MaufaulDa, 5 = others)
T_SEC = T_WAVEFORM
NFFT = 2**10   # (2**10 = MaufaulDa, 2**14 = others)
F_LIMIT = None # (3000 = MaufaulDa, None = others)

import os
import pandas as pd
import numpy as np
import matplotlib.pylab as plt

from tabulate import tabulate
from IPython.display import Markdown, HTML
from tqdm.notebook import tqdm

from typing import List, Tuple
from scipy.signal import find_peaks, butter, lfilter
from tsfel.feature_extraction.features import fundamental_frequency

from zipfile import ZipFile
import sys
sys.path.append('../')
from vibrodiagnostics import mafaulda, selection, discovery, models


def beaglebone_measurement(filename: str, fs: int) -> Tuple[str, pd.DataFrame]:
    g = 9.81
    milivolts = 1800
    resolution = 2**12
    columns = ['x', 'y', 'z']
    ts = pd.read_csv(filename, delimiter='\t', index_col=False, header=None, names=columns)
        
    # Calculate amplitude in m/s^2 Beaglebone Black ADC and ADXL335 resolution (VIN 1.8V, 12bits)
    for dim in columns:
        ts[dim] = ts[dim] * (milivolts / resolution)  # ADC to mV
        ts[dim] = (ts[dim] / 180) * g                 # mV to m/s^2 (180 mV/g)
        ts[dim] -= ts[dim].mean()

    ts['t'] = ts.index * (1 / fs)
    ts.set_index('t', inplace=True)
    return (os.path.basename(filename), ts, fs, ts.columns)  # last is feature columns


def beaglebone_dataset(filenames: List[str], fs: int) -> List[Tuple[str, pd.DataFrame]]:
    dataset = []
    for filename in filenames:
        name, ts, fs, cols = beaglebone_measurement(filename, fs)
        dataset.append((name, ts))
    return dataset


def lowpass_filter(data, cutpoint, fs, order=5):
    b, a = butter(order, cutpoint, fs=fs, btype='lowpass')
    y = lfilter(b, a, data)
    return y


def mafaulda_dataset(
        place=['ax', 'ay', 'az'],
        features_path =  '../../datasets/features_data/',
        mafaulda_path='../../datasets/MAFAULDA.zip',
        rpm=2500,
        lowpass_hz=10000):

    metadata_filename = os.path.join(features_path, selection.MAFAULDA_METADATA)
    faults = {
        'normal': 'normal',
        'imbalance': 'imbalance',
        'horizontal-misalignment': 'misalignment',
        'vertical-misalignment': 'misalignment',
        'overhang-cage_fault': 'cage fault',
        'underhang-cage_fault': 'cage fault',
        'underhang-ball_fault': 'ball fault',
        'overhang-ball_fault': 'ball fault',
        'overhang-outer_race': 'outer race fault',
        'underhang-outer_race': 'outer race fault',
    }

    metadata = pd.read_csv(metadata_filename, index_col='filename')
    metadata.reset_index(inplace=True)
    metadata = models.fault_labeling(metadata, faults)
    files = pd.DataFrame()
    # Worst severity and mid rpm
    for name, group in metadata[metadata['rpm'] >= rpm].groupby(by='fault', observed=False):
        files = pd.concat([
            files,
            group[
                group['severity_level'] == group['severity_level'].max()
            ].sort_values(by='rpm', ascending=True).head(1)
        ])
    ordering = {
        'normal': 0,
        'misalignment': 1,
        'imbalance': 2,
        'cage fault': 3,
        'ball fault': 4,
        'outer race fault': 5,
    }
    source = ZipFile(mafaulda_path)
    dataset = len(files) * [0]
    for index, file in files.iterrows():
        ts = mafaulda.csv_import(source, file['filename'])
        ts = ts[place]
        ts.columns = ts.columns.str.extract(r'(\w)$')[0]
        for axis in ts.columns:
            ts[axis] = lowpass_filter(ts[axis], lowpass_hz, file['fs'])
        pos = ordering[file['fault']]
        dataset[pos] = ((file['fault'] + ' (' + file['filename'] +')', ts))

    return dataset

Load dataset

if DEVICE == 'beaglebone-fan':
    Fs = 2500
    path = '../../inspections/fan/'
    files = [
        '1_still.tsv', '2_still.tsv', '3_still.tsv',
        '1_up.tsv', '2_up.tsv', '3_up.tsv',
        '1_down.tsv', '2_down.tsv', '3_down.tsv'
    ]
    files = [os.path.join(path, name) for name in files]
    DATASET = beaglebone_dataset(files, Fs)

elif DEVICE == 'beaglebone-compressor':
    Fs = 2500
    path = '../../inspections/datacentres/shc3/'
    files = [
        'k3_1.tsv', 'k3_2.tsv', 'k3_3.tsv', 'k3_4.tsv',
        'k5_1.tsv', 'k5_2.tsv', 'k5_3.tsv', 'k5_4.tsv'
    ]
    files = [os.path.join(path, name) for name in files]
    DATASET = beaglebone_dataset(files, Fs)

elif DEVICE == 'beaglebone-pump':
    Fs = 2500
    path = '../../inspections/bvs/'
    files = [
        'bvs_1_hore.tsv', 'bvs_2_hore.tsv' 
        #, 'bvs_3_motor.tsv', 'bvs_4_motor.tsv'
    ]    
    files = [os.path.join(path, name) for name in files]
    DATASET = beaglebone_dataset(files, Fs)

elif DEVICE == 'beaglebone-refrigerator':
    Fs = 2500
    path = '../../inspections/home-refrigerator/'
    files = [
        'ch1.tsv', 'ch2.tsv', 'ch3.tsv', 'ch4.tsv', 'ch5.tsv'
    ]    
    files = [os.path.join(path, name) for name in files]
    DATASET = beaglebone_dataset(files, Fs)

elif DEVICE == 'mafaulda-a':
    Fs = mafaulda.FS_HZ
    DATASET = mafaulda_dataset(place=['ax', 'ay', 'az'])

elif DEVICE == 'mafaulda-b':
    Fs = mafaulda.FS_HZ
    DATASET = mafaulda_dataset(place=['bx', 'by', 'bz'])

for name, ts in DATASET:
    display(Markdown(f'**{name}**'))
    ts.info()
    print()

k3_1.tsv

<class 'pandas.core.frame.DataFrame'>
Index: 51000 entries, 0.0 to 20.3996
Data columns (total 3 columns):
 #   Column  Non-Null Count  Dtype  
---  ------  --------------  -----  
 0   x       51000 non-null  float64
 1   y       51000 non-null  float64
 2   z       51000 non-null  float64
dtypes: float64(3)
memory usage: 1.6 MB

k3_2.tsv

<class 'pandas.core.frame.DataFrame'>
Index: 51000 entries, 0.0 to 20.3996
Data columns (total 3 columns):
 #   Column  Non-Null Count  Dtype  
---  ------  --------------  -----  
 0   x       51000 non-null  float64
 1   y       51000 non-null  float64
 2   z       51000 non-null  float64
dtypes: float64(3)
memory usage: 1.6 MB

k3_3.tsv

<class 'pandas.core.frame.DataFrame'>
Index: 51000 entries, 0.0 to 20.3996
Data columns (total 3 columns):
 #   Column  Non-Null Count  Dtype  
---  ------  --------------  -----  
 0   x       51000 non-null  float64
 1   y       51000 non-null  float64
 2   z       51000 non-null  float64
dtypes: float64(3)
memory usage: 1.6 MB

k3_4.tsv

<class 'pandas.core.frame.DataFrame'>
Index: 51000 entries, 0.0 to 20.3996
Data columns (total 3 columns):
 #   Column  Non-Null Count  Dtype  
---  ------  --------------  -----  
 0   x       51000 non-null  float64
 1   y       51000 non-null  float64
 2   z       51000 non-null  float64
dtypes: float64(3)
memory usage: 1.6 MB

k5_1.tsv

<class 'pandas.core.frame.DataFrame'>
Index: 51000 entries, 0.0 to 20.3996
Data columns (total 3 columns):
 #   Column  Non-Null Count  Dtype  
---  ------  --------------  -----  
 0   x       51000 non-null  float64
 1   y       51000 non-null  float64
 2   z       51000 non-null  float64
dtypes: float64(3)
memory usage: 1.6 MB

k5_2.tsv

<class 'pandas.core.frame.DataFrame'>
Index: 51000 entries, 0.0 to 20.3996
Data columns (total 3 columns):
 #   Column  Non-Null Count  Dtype  
---  ------  --------------  -----  
 0   x       51000 non-null  float64
 1   y       51000 non-null  float64
 2   z       51000 non-null  float64
dtypes: float64(3)
memory usage: 1.6 MB

k5_3.tsv

<class 'pandas.core.frame.DataFrame'>
Index: 51000 entries, 0.0 to 20.3996
Data columns (total 3 columns):
 #   Column  Non-Null Count  Dtype  
---  ------  --------------  -----  
 0   x       51000 non-null  float64
 1   y       51000 non-null  float64
 2   z       51000 non-null  float64
dtypes: float64(3)
memory usage: 1.6 MB

k5_4.tsv

<class 'pandas.core.frame.DataFrame'>
Index: 51000 entries, 0.0 to 20.3996
Data columns (total 3 columns):
 #   Column  Non-Null Count  Dtype  
---  ------  --------------  -----  
 0   x       51000 non-null  float64
 1   y       51000 non-null  float64
 2   z       51000 non-null  float64
dtypes: float64(3)
memory usage: 1.6 MB

for name, ts in DATASET:
    display(Markdown(f'**{name}**'))
    display(tabulate(ts.describe(), headers='keys', tablefmt='html'))
    ts.boxplot(grid=True)
    plt.show()

k3_1.tsv

	x	y	z
count	51000	51000	51000
mean	2.30926e-15	1.14551e-15	3.23924e-16
std	1.37145	2.11749	0.512991
min	-2.46212	-5.23984	-1.22658
25%	-1.31252	-1.59941	-0.436227
50%	0.076596	-0.377948	-0.0530239
75%	1.39386	1.32252	0.47388
max	2.68717	4.69949	1.26424

k3_2.tsv

	x	y	z
count	51000	51000	51000
mean	-1.12739e-15	3.32283e-15	5.67347e-15
std	1.37462	1.90235	0.54732
min	-2.44795	-4.07667	-1.19419
25%	-1.32229	-1.4661	-0.475688
50%	0.0907731	-0.436238	0.00331546
75%	1.43198	1.07262	0.53022
max	2.48579	4.1143	1.36848

k3_3.tsv

	x	y	z
count	51000	51000	51000
mean	9.60069e-16	-1.59607e-15	3.9874e-16
std	1.40137	1.95904	0.55607
min	-2.62483	-4.38867	-1.25534
25%	-1.25967	-1.56255	-0.488938
50%	0.0336402	-0.269236	0.0140165
75%	1.4467	1.26358	0.540921
max	2.76396	4.16155	1.28338

k3_4.tsv

	x	y	z
count	51000	51000	51000
mean	-2.11937e-15	4.49606e-15	5.62652e-15
std	0.627253	0.8348	0.22793
min	-9.79898	-6.27935	-1.67302
25%	-0.0273018	-0.0283453	-0.0444018
50%	0.0205986	-0.0043951	-0.0204516
75%	0.0445488	0.0435053	0.0274488
max	7.75651	9.81518	1.72791

k5_1.tsv

	x	y	z
count	51000	51000	51000
mean	3.58434e-15	-2.53483e-15	6.39767e-16
std	1.49955	2.32511	0.231579
min	-3.22158	-5.17207	-0.665095
25%	-1.09002	-2.01065	-0.162141
50%	0.251196	-0.382032	-0.0184398
75%	1.2571	2.15669	0.149212
max	2.88572	4.95886	1.25092

k5_2.tsv

	x	y	z
count	51000	51000	51000
mean	-9.78738e-16	2.91908e-15	2.61535e-15
std	1.5101	2.32367	0.22692
min	-3.3225	-4.90644	-0.651701
25%	-1.09513	-2.05637	-0.148747
50%	0.270034	-0.451703	-0.0289962
75%	1.29989	2.13492	0.138655
max	2.97641	4.79339	1.19246

k5_3.tsv

	x	y	z
count	51000	51000	51000
mean	-5.92676e-15	-7.54011e-16	1.70753e-15
std	1.49944	2.33441	0.225151
min	-3.29453	-4.94197	-0.632079
25%	-1.09111	-2.044	-0.153075
50%	0.298001	-0.415387	-0.00937392
75%	1.25601	2.12333	0.134327
max	2.76487	4.90156	1.21209

k5_4.tsv

	x	y	z
count	51000	51000	51000
mean	8.0361e-15	8.27908e-15	4.04591e-16
std	1.40058	2.10618	0.258499
min	-4.15798	-5.68461	-1.33782
25%	-1.02051	-1.85257	-0.140307
50%	0.0572508	0.0155409	-0.0445065
75%	1.20686	1.85971	0.123145
max	5.87715	6.93715	3.76357

Statistical tests

• Normality test: Kolmogorov–Smirnov test
• Normality visual test: Quantile-quantile plot on chosen recording
• Stationarity test: Augmented Dickey–Fuller test
• Stationarity visual test: Autocorrelation plot

from statsmodels.tsa.stattools import adfuller
from statsmodels.api import qqplot
from scipy.stats import kstest

normality_tests = []
for name, ts in DATASET:
    for x in ts.columns:
        p_value = kstest(ts[x], 'norm').pvalue
        test = {'name': name, 'axis': x, 'p-value': p_value, 'not-normal': p_value < 0.05}
        normality_tests.append(test)

normality_tests = pd.DataFrame.from_records(normality_tests)
print(normality_tests.value_counts('not-normal'))
normality_tests.describe()

not-normal
True    24
Name: count, dtype: int64

	p-value
count	24.0
mean	0.0
std	0.0
min	0.0
25%	0.0
50%	0.0
75%	0.0
max	0.0

name, ts = DATASET[0]
fig, ax = plt.subplots(1, len(ts.columns), figsize=(10, 4))
for i, x in enumerate(ts.columns):
    qqplot(ts[x], line='45', ax=ax[i], marker='.', alpha=0.5)
    ax[i].set_title(f'Axis: {x}')

plt.tight_layout()
print(name)
plt.show()

k3_1.tsv

stationarity_tests = []
for name, ts in tqdm(DATASET):
    for x in ts.columns:
        result = adfuller(ts[x].loc[T_WAVEFORM:T_WAVEFORM+1])
        p_value = result[1]
        test = {
            'name': name,
            'axis': x,
            'statistic': result[0],
            'p-value': p_value,
            'stationary': p_value < 0.001
        }
        stationarity_tests.append(test)

stationarity_tests = pd.DataFrame.from_records(stationarity_tests)
print(stationarity_tests.value_counts('stationary'))
stationarity_tests['p-value'].describe()

  0%|          | 0/8 [00:00<?, ?it/s]

stationary
True     23
False     1
Name: count, dtype: int64

count    2.400000e+01
mean     7.193267e-04
std      3.523967e-03
min      0.000000e+00
25%      0.000000e+00
50%      1.006159e-29
75%      2.375520e-20
max      1.726384e-02
Name: p-value, dtype: float64

name, ts = DATASET[0]
fig, ax = plt.subplots(1, len(ts.columns), figsize=(10, 4))
for i, x in enumerate(ts.columns):
    ax[i].acorr(ts[x], maxlags=50)
    ax[i].set_title(f'Axis: {x}')

plt.tight_layout()
print(name)
plt.show()

k3_1.tsv

Time domain histogram

for name, ts in DATASET:
    display(Markdown(f'**{name}**'))
    axis = ts.columns
    ax = ts[axis].hist(figsize=(15, 3), grid=True, bins=100, layout=(1, 3), edgecolor='black', linewidth=0.5)
    plt.show()

k3_1.tsv

k3_2.tsv

k3_3.tsv

k3_4.tsv

k5_1.tsv

k5_2.tsv

k5_3.tsv

k5_4.tsv

Time domain waveform

for name, ts in DATASET:
    display(Markdown(f'**{name}**'))
    axis = ts.columns
    
    ax = ts[axis].plot(figsize=(20, 8), grid=True, subplots=True)
    for i, axname in enumerate(axis):
        ax[i].set_xlabel('Time [s]')
        ax[i].set_ylabel(f'Amplitude ({axname}) [m/s^2]')
    plt.show()               # plt.savefig('waveform.png')

k3_1.tsv

k3_2.tsv

k3_3.tsv

k3_4.tsv

k5_1.tsv

k5_2.tsv

k5_3.tsv

k5_4.tsv

Time domain waveform zoom detail

for name, ts in DATASET:
    axis = ts.columns
    display(Markdown(f'**{name}**'))
    ax = (ts[axis].iloc[int(T_WAVEFORM*Fs):int(T_WAVEFORM*Fs)+Fs]
                  .plot(figsize=(20, 10), grid=True, subplots=True))
    
    for i, axname in enumerate(axis):
        ax[i].set_xlabel('Time [s]')
        ax[i].set_ylabel(f'Amplitude ({axname}) [m/s^2]')
        plt.show()      # plt.savefig('waveform_zoom.png')

k3_1.tsv

k3_2.tsv

k3_3.tsv

k3_4.tsv

k5_1.tsv

k5_2.tsv

k5_3.tsv

k5_4.tsv

Time domain waveform zoom - faults side by side

fig, ax = plt.subplots(len(DATASET), 3, figsize=(12, 15), sharex=True)

for idx, df in enumerate(DATASET):
    name, ts = df
    columns = ts.columns
    ax[idx][1].set_title(name)
    ax[idx][0].set_ylabel('Amplitude [m/s^2]')

    for pos, axis in enumerate(columns):
        data = ts[axis].loc[T_WAVEFORM:T_WAVEFORM+0.3]
        ax[idx][pos].plot(data.index, data, linewidth=1, color='darkblue')
        ax[idx][pos].grid()
    

plt.tight_layout()
plt.show()

def spectogram(x, debug=True):
    fig, ax = plt.subplots(figsize=(15, 4))
    cmap = plt.get_cmap('inferno')
    pxx, freqs, t, im = plt.specgram(
        x, NFFT=NFFT, Fs=Fs,
        detrend='mean',
        mode='magnitude', scale='dB',
        cmap=cmap, vmin=-60
    )
    fig.colorbar(im, aspect=20, pad=0.04)
    ax.set_xlabel('Time [s]')
    ax.set_ylabel('Frequency [Hz]')
    mafaulda.resolution_calc(Fs, NFFT)
    return freqs, pxx


def window_idx(t):
    return (Fs * t) // NFFT + 1


def spectrum_slice(freqs, Pxx, t):
    fig, ax = plt.subplots(2, 1, figsize=(20, 8))
    n = window_idx(t)

    dB = 20 * np.log10(Pxx.T[n] / 0.000001)
    ax[0].plot(freqs, dB)      # 1 dB = 1 um/s^2
    ax[0].grid(True)
    ax[0].set_xlabel('Frequency [Hz]')
    ax[0].set_ylabel('Amplitude [dB]')
    
    ax[1].plot(freqs, Pxx.T[n])
    ax[1].grid(True)
    ax[1].set_xlabel('Frequency [Hz]')
    ax[1].set_ylabel('Amplitude [m/s^2]')
    return n


def get_max_frequency(freqs, Pxx, i):
    max_freq = freqs[np.argmax(Pxx.T[i])]
    return max_freq


def get_peaks(freqs, Pxx, i, top=5):
    amplitudes = Pxx.T[i]
    peaks, _ = find_peaks(amplitudes, distance=3)

    fundamental = get_max_frequency(freqs, Pxx, i)
    f_top = freqs[peaks[np.argsort(amplitudes[peaks])]][::-top]
    y_top = np.sort(amplitudes[peaks])[::-top]

    return pd.DataFrame({
        'f': f_top,
        'y': y_top,
        '1x': f_top / fundamental 
    })


def butter_bandpass_filter(data, lowcut, highcut, fs, order=5):
    b, a = butter(order, [lowcut, highcut], fs=fs, btype='band')
    y = lfilter(b, a, data)
    return y


def get_spectrograms(DATASET: List[pd.DataFrame], axis: str) -> list:
    spectrograms = []

    for name, ts in DATASET:
        base_freq = fundamental_frequency(ts[axis], Fs)
        display(Markdown(f'**{name}** *({axis.upper()} axis, Fundamental = {base_freq:.4f} Hz)*'))
        
        freqs, Pxx = spectogram(ts[axis])
        spectrograms.append((name, freqs, Pxx))
        plt.show()          # plt.savefig(f'x_axis_fft_{NFFT}.png')
    
    return spectrograms


def show_spectrogram_detail(spectrograms: list, axis: str, t: float):
    for name, freqs, Pxx in spectrograms:
        display(Markdown(f'**{name}** ({axis.upper()} axis @ {t}s)'))
        i_window = spectrum_slice(freqs, Pxx, t)
        plt.show()           #plt.savefig(f'x_axis_fft_{NFFT}_at_{T_SEC}s.png')


def show_mms_peaks(spectrograms: list, axis: str, t: float):
    for name, freqs, Pxx in spectrograms:
        display(Markdown(f'**{name}** ({axis.upper()} axis @ {t}s)'))
    
        i_window = window_idx(t)
        peaks = discovery.mms_peak_finder(Pxx.T[i_window])
        
        fig, ax = plt.subplots(1, 1, figsize=(15, 3))
        ax.grid(True)
        ax.plot(freqs, Pxx.T[i_window])
        ax.scatter(freqs[peaks], Pxx.T[i_window][peaks], marker='^', color='red')
        ax.set_xlabel('Frequency [Hz]')
        
        plt.show()


def show_harmonic_series(spectrograms: list, axis: str, t: float):
    # https://stackoverflow.com/questions/1982770/changing-the-color-of-an-axis
    for name, freqs, Pxx in spectrograms:
        display(Markdown(f'**{name}** ({axis.upper()} axis @ {t}s)'))
    
        i_window = window_idx(t)
        h_series = discovery.harmonic_series_detection(freqs, Pxx.T[i_window], Fs, NFFT)
    
        # Find best (sum of harmonics' amplitudes in the largest)
        max_harmonic_amp_idx = np.argmax([
            sum([h[1] for h in s]) / len(s)
            for s in h_series
        ])
        best_harmonic_series = pd.DataFrame(
            h_series[max_harmonic_amp_idx],
            columns=['Frequency [Hz]', 'Amplitude [m/s^2]']
        )
        best_harmonic_series.index += 1
        display(tabulate(best_harmonic_series, headers='keys', tablefmt='html'))
    
        # Plot found harmonic series
        fig, ax = plt.subplots(1, 8, figsize=(30, 4))
        for i in range(8):
            s = h_series[i+1]
            if i == max_harmonic_amp_idx:
                ax[i].xaxis.label.set_color('red')
    
            ax[i].plot(freqs, Pxx.T[i_window])
            ax[i].scatter([x[0] for x in s], [x[1] for x in s], marker='^', color='red')
            ax[i].set_xlabel('Frequency [Hz]')
    
        plt.show()

def show_spectra_largest_amplitudes(spectrograms: list, axis: str, t: float):
    for name, freqs, Pxx in spectrograms:
        display(Markdown(f'**{name}** ({axis.upper()} axis @ {t}s)'))

        i_window = window_idx(t)
        x_fundamental = get_max_frequency(freqs, Pxx, i_window)
        peaks = get_peaks(freqs, Pxx, i_window)
        
        display(Markdown(f'- *Fundamental frequency:* {x_fundamental} Hz'))
        display(tabulate(peaks.head(5), headers='keys', tablefmt='html'))


def compare_limited_specrograms(spectrograms: list, axis: str, t: float):
    fig, ax = plt.subplots(len(DATASET), 1, figsize=(20, 20), sharey=True)
    i = 0
    for name, ts in DATASET:
        signal = ts[axis].loc[t:t+NFFT/Fs].to_numpy()
        pxx = np.abs(np.fft.rfft(signal) / len(signal)) 
        freqs = np.fft.fftfreq(len(signal), d=1/Fs)[:len(pxx)]
        #ilast = len(freqs[freqs < F_LIMIT])
        
        ax[i].plot(freqs, pxx)
        ax[i].grid(True)
        ax[i].set_xlabel('Frequency [Hz]')
        ax[i].set_ylabel('Amplitude [m/s^2]')
        ax[i].set_xlim(0, F_LIMIT)
        ax[i].set_ylim(0, 0.4)
        ax[i].set_title(name)
        i += 1


def spectrogram_energy_left_cumulative(spectrograms: list, axis: str, t: float):
    fig, ax = plt.subplots(len(DATASET), 1, figsize=(20, 20), sharey=True)
    i = 0
    for name, ts in DATASET:
        signal = ts[axis].loc[t:t+NFFT/Fs].to_numpy()
        pxx = np.abs(np.fft.rfft(signal) / len(signal)) 
        freqs = np.fft.fftfreq(len(signal), d=1/Fs)[:len(pxx)]
        
        ax[i].plot(freqs, np.cumsum(pxx) / np.sum(pxx))
        ax[i].grid(True)
        ax[i].set_xlabel('Frequency [Hz]')
        ax[i].set_ylabel('Cumulative energy [%]')
        #ax[i].set_xlim(0, 10000)
        ax[i].set_title(name)
        i += 1

Compare mafaulda faults

compare_limited_specrograms(DATASET, 'x', T_SEC)
plt.tight_layout()
plt.show()

compare_limited_specrograms(DATASET, 'y', T_SEC)
plt.tight_layout()
plt.show()

compare_limited_specrograms(DATASET, 'z', T_SEC)
plt.tight_layout()
plt.show()

Compare cumulative sums

spectrogram_energy_left_cumulative(DATASET, 'x', T_SEC)
plt.tight_layout()
plt.show()

spectrogram_energy_left_cumulative(DATASET, 'y', T_SEC)
plt.tight_layout()
plt.show()

spectrogram_energy_left_cumulative(DATASET, 'z', T_SEC)
plt.tight_layout()
plt.show()

Spectrogram in X axis

x_spectra = get_spectrograms(DATASET, 'x')

k3_1.tsv (X axis, Fundamental = 47.7451 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz

k3_2.tsv (X axis, Fundamental = 47.7451 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz

k3_3.tsv (X axis, Fundamental = 47.7451 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz

k3_4.tsv (X axis, Fundamental = 47.7941 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz

k5_1.tsv (X axis, Fundamental = 48.1373 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz

k5_2.tsv (X axis, Fundamental = 48.1373 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz

k5_3.tsv (X axis, Fundamental = 48.1373 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz

k5_4.tsv (X axis, Fundamental = 48.1373 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz

Spectrogram detail in X axis

show_spectrogram_detail(x_spectra, 'x', T_SEC)

k3_1.tsv (X axis @ 5s)

k3_2.tsv (X axis @ 5s)

k3_3.tsv (X axis @ 5s)

k3_4.tsv (X axis @ 5s)

k5_1.tsv (X axis @ 5s)

k5_2.tsv (X axis @ 5s)

k5_3.tsv (X axis @ 5s)

k5_4.tsv (X axis @ 5s)

Peaks in frequency spectrum in X axis

• MMS peak finder algorithm

show_mms_peaks(x_spectra, 'x', T_SEC)

k3_1.tsv (X axis @ 5s)

k3_2.tsv (X axis @ 5s)

k3_3.tsv (X axis @ 5s)

k3_4.tsv (X axis @ 5s)

k5_1.tsv (X axis @ 5s)

k5_2.tsv (X axis @ 5s)

k5_3.tsv (X axis @ 5s)

k5_4.tsv (X axis @ 5s)

Harmonic series detection in X axis

# show_harmonic_series(x_spectra, 'x', T_SEC)

show_spectra_largest_amplitudes(x_spectra, 'x', T_SEC)

k3_1.tsv (X axis @ 5s)

• Fundamental frequency: 48.828125 Hz

	f	y	1x
0	48.8281	0.685933	1
1	285.645	0.180302	5.85
2	524.902	0.0620602	10.75
3	1052.25	0.0087056	21.55
4	717.773	0.00689525	14.7

k3_2.tsv (X axis @ 5s)

• Fundamental frequency: 48.828125 Hz

	f	y	1x
0	48.8281	0.703858	1
1	285.645	0.136529	5.85
2	524.902	0.0330748	10.75
3	954.59	0.00578816	19.55
4	14.6484	0.00447683	0.3

k3_3.tsv (X axis @ 5s)

• Fundamental frequency: 48.828125 Hz

	f	y	1x
0	48.8281	0.656088	1
1	285.645	0.149954	5.85
2	144.043	0.0424332	2.95
3	1164.55	0.0068831	23.85
4	639.648	0.00539692	13.1

k3_4.tsv (X axis @ 5s)

• Fundamental frequency: 4.8828125 Hz

	f	y	1x
0	4.88281	0.0768066	1
1	153.809	0.00184415	31.5
2	297.852	0.00152948	61
3	83.0078	0.00119425	17
4	1015.62	0.0010988	208

k5_1.tsv (X axis @ 5s)

• Fundamental frequency: 48.828125 Hz

	f	y	1x
0	48.8281	0.921647	1
1	288.086	0.102319	5.9
2	434.57	0.0270554	8.9
3	722.656	0.0109019	14.8
4	400.391	0.00575804	8.2

k5_2.tsv (X axis @ 5s)

• Fundamental frequency: 48.828125 Hz

	f	y	1x
0	48.8281	0.92099	1
1	288.086	0.10624	5.9
2	480.957	0.0224154	9.85
3	14.6484	0.00968361	0.3
4	627.441	0.00550303	12.85

k5_3.tsv (X axis @ 5s)

• Fundamental frequency: 48.828125 Hz

	f	y	1x
0	48.8281	0.936292	1
1	288.086	0.106108	5.9
2	241.699	0.0264785	4.95
3	295.41	0.0106427	6.05
4	1103.52	0.00580459	22.6

k5_4.tsv (X axis @ 5s)

• Fundamental frequency: 48.828125 Hz

	f	y	1x
0	48.8281	0.898246	1
1	288.086	0.0986729	5.9
2	480.957	0.027933	9.85
3	295.41	0.00986706	6.05
4	983.887	0.00442585	20.15

Spectrogram in Y axis

y_spectra = get_spectrograms(DATASET, 'y')

k3_1.tsv (Y axis, Fundamental = 47.7451 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz

k3_2.tsv (Y axis, Fundamental = 47.7451 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz

k3_3.tsv (Y axis, Fundamental = 47.7451 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz

k3_4.tsv (Y axis, Fundamental = 47.3039 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz

k5_1.tsv (Y axis, Fundamental = 48.1373 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz

k5_2.tsv (Y axis, Fundamental = 48.1373 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz

k5_3.tsv (Y axis, Fundamental = 48.1373 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz

k5_4.tsv (Y axis, Fundamental = 48.1373 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz

Spectrogram detail in Y axis

show_spectrogram_detail(y_spectra, 'y', T_SEC)

k3_1.tsv (Y axis @ 5s)

k3_2.tsv (Y axis @ 5s)

k3_3.tsv (Y axis @ 5s)

k3_4.tsv (Y axis @ 5s)

k5_1.tsv (Y axis @ 5s)

k5_2.tsv (Y axis @ 5s)

k5_3.tsv (Y axis @ 5s)

k5_4.tsv (Y axis @ 5s)

Peaks in frequency spectrum in Y axis

show_mms_peaks(y_spectra, 'y', T_SEC)

k3_1.tsv (Y axis @ 5s)

k3_2.tsv (Y axis @ 5s)

k3_3.tsv (Y axis @ 5s)

k3_4.tsv (Y axis @ 5s)

k5_1.tsv (Y axis @ 5s)

k5_2.tsv (Y axis @ 5s)

k5_3.tsv (Y axis @ 5s)

k5_4.tsv (Y axis @ 5s)

Harmonic series detection in Y axis

# show_harmonic_series(y_spectra, 'y', T_SEC)

show_spectra_largest_amplitudes(y_spectra, 'y', T_SEC)

k3_1.tsv (Y axis @ 5s)

• Fundamental frequency: 48.828125 Hz

	f	y	1x
0	48.8281	1.20267	1
1	478.516	0.198349	9.8
2	285.645	0.103564	5.85
3	573.73	0.0188074	11.75
4	622.559	0.0106902	12.75

k3_2.tsv (Y axis @ 5s)

• Fundamental frequency: 48.828125 Hz

	f	y	1x
0	48.8281	1.03286	1
1	144.043	0.138609	2.95
2	285.645	0.0364966	5.85
3	541.992	0.0092398	11.1
4	1240.23	0.00614032	25.4

k3_3.tsv (Y axis @ 5s)

• Fundamental frequency: 48.828125 Hz

	f	y	1x
0	48.8281	1.00388	1
1	478.516	0.12413	9.8
2	429.688	0.0633835	8.8
3	1164.55	0.0102918	23.85
4	1115.72	0.00770109	22.85

k3_4.tsv (Y axis @ 5s)

• Fundamental frequency: 9.765625 Hz

	f	y	1x
0	9.76562	0.0299383	1
1	21.9727	0.00412541	2.25
2	83.0078	0.00197565	8.5
3	244.141	0.00151702	25
4	227.051	0.0013613	23.25

k5_1.tsv (Y axis @ 5s)

• Fundamental frequency: 48.828125 Hz

	f	y	1x
0	48.8281	1.38965	1
1	385.742	0.147958	7.9
2	529.785	0.0615459	10.85
3	400.391	0.0299462	8.2
4	26.8555	0.0151213	0.55

k5_2.tsv (Y axis @ 5s)

• Fundamental frequency: 48.828125 Hz

	f	y	1x
0	48.8281	1.40826	1
1	434.57	0.179166	8.9
2	419.922	0.0564781	8.6
3	405.273	0.0231326	8.3
4	673.828	0.0110141	13.8

k5_3.tsv (Y axis @ 5s)

• Fundamental frequency: 48.828125 Hz

	f	y	1x
0	48.8281	1.42221	1
1	480.957	0.15385	9.85
2	578.613	0.0547047	11.85
3	402.832	0.0217529	8.25
4	1149.9	0.0119277	23.55

k5_4.tsv (Y axis @ 5s)

• Fundamental frequency: 48.828125 Hz

	f	y	1x
0	48.8281	1.3829	1
1	480.957	0.143644	9.85
2	578.613	0.0521242	11.85
3	983.887	0.0168376	20.15
4	14.6484	0.0122953	0.3

Spectrogram in Z axis

z_spectra = get_spectrograms(DATASET, 'z')

k3_1.tsv (Z axis, Fundamental = 47.7451 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz

k3_2.tsv (Z axis, Fundamental = 47.7451 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz

k3_3.tsv (Z axis, Fundamental = 47.7451 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz

k3_4.tsv (Z axis, Fundamental = 47.7451 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz

k5_1.tsv (Z axis, Fundamental = 0.0490 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz

k5_2.tsv (Z axis, Fundamental = 48.1373 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz

k5_3.tsv (Z axis, Fundamental = 48.1373 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz

k5_4.tsv (Z axis, Fundamental = 48.1373 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz

Spectrogram detail in Z axis

show_spectrogram_detail(z_spectra, 'z', T_SEC)

k3_1.tsv (Z axis @ 5s)

k3_2.tsv (Z axis @ 5s)

k3_3.tsv (Z axis @ 5s)

k3_4.tsv (Z axis @ 5s)

k5_1.tsv (Z axis @ 5s)

k5_2.tsv (Z axis @ 5s)

k5_3.tsv (Z axis @ 5s)

k5_4.tsv (Z axis @ 5s)

Peaks in frequency spectrum in Z axis

show_mms_peaks(z_spectra, 'z', T_SEC)

k3_1.tsv (Z axis @ 5s)

k3_2.tsv (Z axis @ 5s)

k3_3.tsv (Z axis @ 5s)

k3_4.tsv (Z axis @ 5s)

k5_1.tsv (Z axis @ 5s)

k5_2.tsv (Z axis @ 5s)

k5_3.tsv (Z axis @ 5s)

k5_4.tsv (Z axis @ 5s)

Harmonic series detection in Z axis

# show_harmonic_series(z_spectra, 'z', T_SEC)

show_spectra_largest_amplitudes(z_spectra, 'z', T_SEC)

k3_1.tsv (Z axis @ 5s)

• Fundamental frequency: 48.828125 Hz

	f	y	1x
0	48.8281	0.280503	1
1	444.336	0.0312579	9.1
2	383.301	0.0217816	7.85
3	507.812	0.0143637	10.4
4	859.375	0.0133318	17.6

k3_2.tsv (Z axis @ 5s)

• Fundamental frequency: 48.828125 Hz

	f	y	1x
0	48.8281	0.32505	1
1	429.688	0.0253058	8.8
2	144.043	0.0191175	2.95
3	954.59	0.0163078	19.55
4	1049.8	0.0117508	21.5

k3_3.tsv (Z axis @ 5s)

• Fundamental frequency: 48.828125 Hz

	f	y	1x
0	48.8281	0.323556	1
1	334.473	0.0279629	6.85
2	668.945	0.0219133	13.7
3	456.543	0.0152852	9.35
4	922.852	0.010976	18.9

k3_4.tsv (Z axis @ 5s)

• Fundamental frequency: 4.8828125 Hz

	f	y	1x
0	4.88281	0.0167918	1
1	61.0352	0.00451767	12.5
2	90.332	0.00312322	18.5
3	324.707	0.00191779	66.5
4	424.805	0.00173085	87

k5_1.tsv (Z axis @ 5s)

• Fundamental frequency: 48.828125 Hz

	f	y	1x
0	48.8281	0.0713598	1
1	434.57	0.0264955	8.9
2	192.871	0.0225645	3.95
3	34.1797	0.0192975	0.7
4	1010.74	0.0169037	20.7

k5_2.tsv (Z axis @ 5s)

• Fundamental frequency: 48.828125 Hz

	f	y	1x
0	48.8281	0.0682579	1
1	192.871	0.0260875	3.95
2	336.914	0.0217514	6.9
3	144.043	0.0200225	2.95
4	31.7383	0.0125301	0.65

k5_3.tsv (Z axis @ 5s)

• Fundamental frequency: 48.828125 Hz

	f	y	1x
0	48.8281	0.0773615	1
1	144.043	0.0257784	2.95
2	7.32422	0.0197851	0.15
3	288.086	0.0169507	5.9
4	422.363	0.0140544	8.65

k5_4.tsv (Z axis @ 5s)

• Fundamental frequency: 48.828125 Hz

	f	y	1x
0	48.8281	0.0907712	1
1	144.043	0.0249169	2.95
2	625	0.0199538	12.8
3	964.355	0.0172125	19.75
4	14.6484	0.0139624	0.3

Histogram

axis = ['x', 'y', 'z']
for name, ts in DATASET:
    display(Markdown(f'**{name}**'))
    ts[axis].hist(figsize=(10, 5), grid=True, bins=50)
    plt.show()

k3_1.tsv

k3_2.tsv

k3_3.tsv

k3_4.tsv

k5_1.tsv

k5_2.tsv

k5_3.tsv

k5_4.tsv

axis = ['x', 'y', 'z']
for name, ts in DATASET:
    display(Markdown(f'**{name}**'))
    ts[axis].boxplot(figsize=(10, 5))
    plt.show()

k3_1.tsv

k3_2.tsv

k3_3.tsv

k3_4.tsv

k5_1.tsv

k5_2.tsv

k5_3.tsv

k5_4.tsv

Orbitals of all cross sections

for name, ts in DATASET:
    display(Markdown(f'**{name}**'))
    fig, ax = plt.subplots(1, 3, figsize=(20, 4))

    for i, col in enumerate([('x', 'y'), ('x', 'z'), ('y', 'z')]):
        ax[i].scatter(ts[col[0]], ts[col[1]], s=1)
        ax[i].grid(True)
        ax[i].set_xlabel(col[0].upper())
        ax[i].set_ylabel(col[1].upper())
        ax[i].grid(True)
    plt.show()       # plt.savefig('orbitals.png')

k3_1.tsv

k3_2.tsv

k3_3.tsv

k3_4.tsv

k5_1.tsv

k5_2.tsv

k5_3.tsv

k5_4.tsv

Orbitals of 1x harmonic frequency

x_spectra_by_name = {spec[0]: spec for spec in x_spectra}
y_spectra_by_name = {spec[0]: spec for spec in y_spectra}
z_spectra_by_name = {spec[0]: spec for spec in z_spectra}
t = 5
space = 5

for name, ts in DATASET:
    display(Markdown(f'**{name}**'))
    fig, ax = plt.subplots(1, 3, figsize=(20, 4))

    name, freqs, Pxx = x_spectra_by_name[name]
    x_fundamental = get_max_frequency(freqs, Pxx, window_idx(t))
    name, freqs, Pxx = y_spectra_by_name[name]
    y_fundamental = get_max_frequency(freqs, Pxx, window_idx(t))

    name, freqs, Pxx = z_spectra_by_name[name]
    z_fundamental = get_max_frequency(freqs, Pxx, window_idx(t))

    try:
        ts['x_1x'] = butter_bandpass_filter(ts['x'], x_fundamental - space, x_fundamental + space, Fs)
        ts['y_1x'] = butter_bandpass_filter(ts['y'], y_fundamental - space, y_fundamental + space, Fs)
        ts['z_1x'] = butter_bandpass_filter(ts['z'], z_fundamental - space, z_fundamental + space, Fs)
    except ValueError:
        continue
    
    for i, col in enumerate([('x_1x', 'y_1x'), ('x_1x', 'z_1x'), ('y_1x', 'z_1x')]):
        ax[i].scatter(ts[col[0]], ts[col[1]], s=1)
        ax[i].grid(True)
        ax[i].set_xlabel(col[0].upper())
        ax[i].set_ylabel(col[1].upper())
        ax[i].grid(True)
    
    plt.show()       # plt.savefig('orbitals_1x.png')

k3_1.tsv

k3_2.tsv

k3_3.tsv

k3_4.tsv

k5_1.tsv

k5_2.tsv

k5_3.tsv

k5_4.tsv

t = 5
space = 8

for name, ts in DATASET:
    display(Markdown(f'**{name}**'))

    name, freqs, Pxx = x_spectra_by_name[name]
    x_fundamental = get_max_frequency(freqs, Pxx, window_idx(t))
    name, freqs, Pxx = y_spectra_by_name[name]
    y_fundamental = get_max_frequency(freqs, Pxx, window_idx(t))

    name, freqs, Pxx = z_spectra_by_name[name]
    z_fundamental = get_max_frequency(freqs, Pxx, window_idx(t))

    try:
        x = butter_bandpass_filter(ts['x'], x_fundamental - space, x_fundamental + space, Fs)
        y = butter_bandpass_filter(ts['y'], y_fundamental - space, y_fundamental + space, Fs)
        z = butter_bandpass_filter(ts['z'], z_fundamental - space, z_fundamental + space, Fs)
    except ValueError:
        continue

    ax = plt.figure().add_subplot(projection='3d')
    ax.scatter(x, y, z, zdir='z', s=1, color='navy')
    ax.set_xlabel('x')
    ax.set_ylabel('y')
    ax.set_zlabel('z')
    ax.set_xlim(-2, 2)
    ax.set_ylim(-2, 2)
    ax.set_zlim(-2, 2)
    ax.zaxis.labelpad = -0.7
    plt.show()

k3_1.tsv

k3_2.tsv

k3_3.tsv

k3_4.tsv

k5_1.tsv

k5_2.tsv

k5_3.tsv

k5_4.tsv